function [W, L] = compute_geo_laplacian(X, knn, wflag, sigma)
    if nargin < 3
        wflag = 'heat';
    end
    
    if nargin < 4
        sigma = 1;
    end
    
    %%% compute adjacency matrix %%%
    np = size(X, 1);
    [nn, dt] = knnsearch(X, X, 'k', knn + 1);
    nn = nn';
    dt = dt';

    Wt = ones(size(nn));
    i = repmat(1: np, knn + 1, 1);
    Wt = sparse(i(:), double(nn(:)), Wt(:), np, np);
    Wt = max(Wt, Wt'); % make symmetric %

    %%% compute weight matrix %%%
    switch wflag
        case 'simple'
            W = Wt;
        case 'heat'
            Gt = sparse(i(:), double(nn(:)), dt(:), np, np);
            Gt = max(Gt, Gt') / max(max(Gt));
            W = Gt .^ 2;
            W(W ~= 0) = exp(-W(W ~= 0) / (2 * sigma ^ 2));
    end

    if nargout > 1
        %%% compute laplacian & eigen-decomp matrix %%%
        Dh = diag(sum(W, 2) .^ -0.5);
        Dn = Dh .^ 2;
        Wg = Dn * W * Dn;
        Dg = diag(sum(Wg, 2) .^ -1);
        tmp = Dg * Wg;
        tmpt = 1: (size(tmp, 1) + 1): numel(tmp);
        tmp(tmpt) = tmp(tmpt) - 1;
        L = tmp;
    end
end